#################################################################################################################
# Data for the replication of RADIO PUBLIC SERVICE ANNOUNCEMENTS AND VOTER PARTICIPATION AMONG NATIVE AMERICANS: 
#   EVIDENCE FROM TWO FIELD EXPERIMENTS
# Eline A. de Rooij and Donald P. Green
# Political Behavior (POBE)
# 
# FIGURE 1 MANUSCRIPT: Radio Bayes
# Code written by Alex Coppock and updated by Don Green 7-8-16
#################################################################################################################

rm(list=ls())

library(ggplot2)
library(dplyr)

# Helper Functions

bayes_updater <- function(est_1, se_1, est_2, se_2){
  est_pooled <- weighted.mean(c(est_1, est_2), w = 1/(c(se_1, se_2)^2))
  se_pooled <- sqrt(1/((1/se_1^2) + (1/se_2^2)))
  return(c(est_pooled, se_pooled))
}

radio_learning <- function(prior_mean, prior_sd){
  update_1 <- bayes_updater(est_1 = prior_mean, se_1 = prior_sd,
                            est_2 = estimates[1], se_2 = ses[1])
  
  update_2 <- bayes_updater(est_1 = update_1[1], se_1 = update_1[2],
                            est_2 = estimates[2], se_2 = ses[2])
  
  return(c(prior_mean, prior_sd, update_1, update_2))
}


multiplot <- function(..., plotlist=NULL, file, cols=1, layout=NULL) {
  require(grid)
  
  # Make a list from the ... arguments and plotlist
  plots <- c(list(...), plotlist)
  
  numPlots = length(plots)
  
  # If layout is NULL, then use 'cols' to determine layout
  if (is.null(layout)) {
    # Make the panel
    # ncol: Number of columns of plots
    # nrow: Number of rows needed, calculated from # of cols
    layout <- matrix(seq(1, cols * ceiling(numPlots/cols)),
                     ncol = cols, nrow = ceiling(numPlots/cols))
  }
  
  if (numPlots==1) {
    print(plots[[1]])
    
  } else {
    # Set up the page
    grid.newpage()
    pushViewport(viewport(layout = grid.layout(nrow(layout), ncol(layout))))
    
    # Make each plot, in the correct location
    for (i in 1:numPlots) {
      # Get the i,j matrix positions of the regions that contain this subplot
      matchidx <- as.data.frame(which(layout == i, arr.ind = TRUE))
      
      print(plots[[i]], vp = viewport(layout.pos.row = matchidx$row,
                                      layout.pos.col = matchidx$col))
    }
  }
}


# Data from experiments: results updated 7-8-16

estimates <- c(0.01681381, 0.01986812)
ses <- c(0.02040339, 0.02521854)


scenario_1 <- radio_learning(prior_mean = 0, prior_sd = 0.1)
scenario_2 <- radio_learning(prior_mean = 0.01, prior_sd = 0.02)
scenario_3 <- radio_learning(prior_mean = 0.02, prior_sd = 0.02)

x <- seq(-0.15, .15, 0.0001)

df <- rbind(cbind(x, y = dnorm(x,mean = scenario_1[1], sd = scenario_1[2]), time = "Prior", scenario = "Scenario 1"),
            cbind(x, y = dnorm(x,mean = scenario_1[3], sd = scenario_1[4]), time = "Update 1", scenario = "Scenario 1"),
            cbind(x, y = dnorm(x,mean = scenario_1[5], sd = scenario_1[6]), time = "Update 2", scenario = "Scenario 1"),
            cbind(x, y = dnorm(x,mean = scenario_2[1], sd = scenario_2[2]), time = "Prior", scenario = "Scenario 2"),
            cbind(x, y = dnorm(x,mean = scenario_2[3], sd = scenario_2[4]), time = "Update 1", scenario = "Scenario 2"),
            cbind(x, y = dnorm(x,mean = scenario_2[5], sd = scenario_2[6]), time = "Update 2", scenario = "Scenario 2"),
            cbind(x, y = dnorm(x,mean = scenario_3[1], sd = scenario_3[2]), time = "Prior", scenario = "Scenario 3"),
            cbind(x, y = dnorm(x,mean = scenario_3[3], sd = scenario_3[4]), time = "Update 1", scenario = "Scenario 3"),
            cbind(x, y = dnorm(x,mean = scenario_3[5], sd = scenario_3[6]), time = "Update 2", scenario = "Scenario 3"))

df_positive <- 
  data.frame(means = c(scenario_1[1], scenario_1[3], scenario_1[5], 
                       scenario_2[1], scenario_2[3], scenario_2[5], 
                       scenario_3[1], scenario_3[3], scenario_3[5]),
             sds = c(scenario_1[2], scenario_1[4], scenario_1[6],   
                       scenario_2[2], scenario_2[4], scenario_2[6], 
                       scenario_3[2], scenario_3[4], scenario_3[6]),
             time = rep(c("Prior", "Update 1", "Update 2"), 3),
             scenario = rep(c("Scenario 1", "Scenario 2", "Scenario 3"), each = 3))

df_joined <- 
  left_join(data.frame(df), df_positive) %>%
  mutate(prop_positive = pnorm(0, mean = means , sd = sds, lower.tail = FALSE),
         facet = paste0(time, ": N(", sprintf("%.3f", means), ",",sprintf("%.3f", sds), ")",
                        "\n pr(ATE > 0) = ", sprintf("%.3f", prop_positive)),
         x=as.numeric(as.character(x)),
         y=as.numeric(as.character(y)))


scenario_1_df <- filter(df_joined, scenario=="Scenario 1")
scenario_2_df <- filter(df_joined, scenario=="Scenario 2")
scenario_3_df <- filter(df_joined, scenario=="Scenario 3")


g1 <- 
ggplot(scenario_1_df, aes(x=x, y=y)) + 
  geom_line() + geom_vline(xintercept = 0, linetype="dashed") + 
  facet_grid(scenario~facet) + 
  theme_bw() + 
  geom_ribbon(data=subset(scenario_1_df, x > 0 ),aes(ymax=y),ymin=0,
              fill="lightgray",colour=NA,alpha=0.5) + 
  theme(axis.text.y=element_blank(),axis.ticks.y=element_blank(),
        axis.title.x=element_blank(),
        axis.title.y=element_blank(),legend.position="none",
        strip.background = element_blank()) + coord_cartesian(xlim=c(-.11, .11))

g2 <- 
  ggplot(scenario_2_df, aes(x=x, y=y)) + 
  geom_line() + geom_vline(xintercept = 0, linetype="dashed") + 
  facet_grid(scenario~facet) + 
  theme_bw() + 
  geom_ribbon(data=subset(scenario_2_df, x > 0 ),aes(ymax=y),ymin=0,
              fill="lightgray",colour=NA,alpha=0.5) + 
  theme(axis.text.y=element_blank(),axis.ticks.y=element_blank(),
        axis.title.x=element_blank(),
        axis.title.y=element_blank(),legend.position="none",
        strip.background = element_blank()) + coord_cartesian(xlim=c(-.11, .11))

g3 <- 
  ggplot(scenario_3_df, aes(x=x, y=y)) + 
  geom_line() + geom_vline(xintercept = 0, linetype="dashed") + 
  facet_grid(scenario~facet) + 
  theme_bw() + 
  geom_ribbon(data=subset(scenario_3_df, x > 0 ),aes(ymax=y),ymin=0,
              fill="lightgray",colour=NA,alpha=0.5) + 
  theme(axis.text.y=element_blank(),axis.ticks.y=element_blank(),
        axis.title.x=element_blank(),
        axis.title.y=element_blank(),legend.position="none",
        strip.background = element_blank()) + coord_cartesian(xlim=c(-.11, .11))



png(file = "Bayesian_Learning_For_Radios.png", width = 900, height = 600)
multiplot(g1, g2, g3)
dev.off()



